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§ ; ABSTRACT 

in ■ 

, We use a multiannulus planetesimal accretion code to investigate the growth of 

I icy planets in the outer regions of a planetesimal disk. In a quiescent minimum mass 

"^H I solar nebula, icy planets grow to sizes of 1000-3000 km on a timescale tp « 15 — 

Dh' 20 Myr (a/30 AU)^, where a is the distance from the central star. Planets form faster 

O ' in more massive nebulae. Newly-formed planets stir up leftover planetesimals along 

' their orbits and produce a collisional cascade where icy planetesimals are slowly ground 

Ci ■ to dust. 

J> . The dusty debris of planet formation has physical characteristics similar to those 

^ , observed in (3 Pic, HR 4796A, and other debris disks. The computed dust masses 

^ : are Md{r < 1 mm) - 10^6 g (Mq/Mmmsn) and Mrf(l mm < r < 1 m) ~ 10^^ g 

■ " " ' (Mq/Mmmsn), where r is the radius of a particle, Mq is the initial mass in solids, and 

Mm MSN is the mass in solids of a minimum mass solar nebula at 30-150 AU. The lu- 
minosity of the dusty disk relative to the stellar luminosity is Ld/Lq ~ Lmax{t/to)~^, 
where Lmax ~ W~^{Mq/Mmmsn), to « 10 Myr to 1 Gyr, and m 1-2. Our calcu- 
lations produce bright rings and dark gaps with sizes Aa/a ~ 0.1. Bright rings occur 
where 1000 km and larger planets have recently formed. Dark gaps are regions where 
planets have cleared out dust or shadows where planets have yet to form. 

Planets can also grow in a planetesimal disk perturbed by the close passage of a star. 
Stellar flybys initiate collisional cascades, which produce copious amounts of dust. The 
dust luminosity following a modest perturbation is 3-4 times larger than the maximum 
dust luminosity of a quiescent planet-forming disk. In 10 Myr or less, large perturbations 
remove almost all of the planetesimals from a disk. After a modest flyby, collisional 
damping reduces planetesimal velocities and allows planets to grow from the remaining 
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planetesimals. Planet formation timescales are then 2-4 times longer than timescales 
for undisturbed disks; dust luminosities are 2-4 times smaller. 

Subject headings: planetary systems - solar system: formation - stars: formation - 
circumstellar matter 



1. INTRODUCTION 

Observations demonstrate that most stars are born with circumstellar disks of gas and dust. 
The disks have masses and sizes similar to the protosolar nebula that formed our solar system 
(Beckwith 1999; Lada 1999). As they age, stars lose their circumstellar disks. In nearby young 
star clusters, 1 Myr old stars often have massive optically thick disks; 10 Myr old stars rarely have 
opaque disks (Haisch, Lada, & Lada 2001). However, older stars often have a 'debris disk' of dusty 
material with a size comparable to the radius of the Kuiper Belt of our solar system (Habing et 
al. 2001; Song et al. 2000; Spangler et al. 2001; Luu & Jewitt 2002). Debris disk masses are 
at least a few lunar masses and may exceed the mass, ~ 0.1 Mq, of the Kuiper Belt (Backman 
& Paresce 1993; Artymowicz 1997; Lagrange et al. 2000; Luu & Jewitt 2002; Wyatt, Dent, & 
Greaves 2003; Greaves & Wyatt 2003). 

Recent imaging and photometric data suggest that debris disks often have internal structure. 
Several disks have asymmetries and warps in their surface brightness distributions (Lagrange et al. 
2000; Augereau et al. 2001). Distinct rings of dust surround (3 Pic (Kalas et al. 2000; Wahhaj 
et al. 2003), e Eri (Greaves et al. 1998), HR 4796A (Jayawardhana ct al. 1998; Koerner et al. 
1998; Augereau et al. 1999; Schneider et al. 1999; Greaves, et al. 2000), and Vega (Wilner et al. 
2002; Koerner, Sargent, & Ostroff 2001). HD 141569 has a dark gap in its disk (Weinberger et al. 
1999). Photometric eclipses in KH 15D suggest at least one discrete, opaque object in orbit around 
the central pre-main sequence star (Herbst et al. 2002). The nearby A- type star Fomalhaut also 
has a distinct clump in its debris disk (Wyatt & Dent 2002; Holland et al. 2003). 

Producing internal structure in a dusty debris disk requires a two (or more) step process. 
Because radiation removes dust from the disk on timescales shorter than the stellar age, some 
process must replenish the dust. Collisions between large bodies form copious amounts of dust if 
the collision velocity is ~ 100-300 m s~^ (e.g. Tanaka et al. 1996; Kenyon et al. 1999; Krivov 
et al. 2000). Small embedded planets (Kenyon et al. 1999) and close encounters with field stars 
(Ida et al. 2000; Kalas et al. 2000) can excite large collision velocities and initiate a 'coUisional 
cascade,' where 1-10 km planetesimals are slowly ground into fine dust grains. The reservoir of 
mass needed to maintain a collisional cascade for the 100-500 Myr lifetimes of most known debris 
disk systems is ~ 10-100 (Habing et al. 2001; Spangler et al. 2001; Wyatt, Dent, & Greaves 
2003) . This mass is ~ 10% to 100% of the amount of solid material in the 'minimum mass solar 
nebula', the minimum amount of material needed to form the planets in our solar system scaled to 
solar abundances (Weidenschilling 1977a; Hayashi 1981). 
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During the collisional cascade, dynamical processes impose structure on the dust. The gravity 
of a passing star can produce ring-hke density concentrations in a particle disk (Kalas et al. 2000). 
Radiation pressure ejects small dust grains and forms dusty rings at heliocentric distances of 80-100 
AU (Takeuchi & Artymowicz 2001). Small planets stir up leftover planetesimals along their orbits 
and produce bright rings of dust (Kenyon & Bromley 2002b). Larger planets or shepherd moons 
can produce dark gaps, partial rings, or warps in the disk (Wilner et al. 2002; Ozernoy et al. 
2000; Wyatt et al. 1999; Wyatt 2003). 

To investigate formation mechanisms for debris disks, we consider dusty disks as plausible 
remnants of recent planet formation. In the planetesimal theory, planets grow from mergers of 
smaller objects embedded in the disk (Safronov 1969). Small dust grains in the disk grow to mm 
sizes and settle into a thin, dusty layer at the midplane of the disk. Collisions between grains in the 
midplane may form successively larger grains which grow into 1 km 'planetesimals' (Weidenschilling 
1980; Weidenschilling & Cuzzi 1993; Kornet, Stepinski, & Rozyczka 2001). Planetesimals may 
form directly through gravitational instabilities (Goldreich &; Ward 1973; Youdin &; Shu 2003). In 
either case, slowly-moving planetesimals collide and merge to form larger and larger bodies which 
eventually become planets (Wetherill & Stewart 1993; Weidenschilling et al. 1997; Kenyon &: Luu 
1999). The gravity of 1000 km or larger planets stirs up leftover planetesimals to large velocities 
(Kenyon &; Bromley 2001). Collisions between rapidly moving planetesimals produce observable 
amounts of dust, which often lie in ring-like structures (Kenyon et al. 1999; Kenyon & Bromley 
2002b). 

In this paper, we continue our study of the bright rings and dark gaps formed during the 
growth of icy planets in the outer regions of a planetesimal disk. Kenyon &; Bromley (2002b) show 
that icy planets reach sizes of 1000-3000 km on timescales of ~ 10-20 Myr at 30 AU and 500 Myr to 
1 Gyr at 100 AU. The dusty debris of icy planet formation often lies in bright rings with lifetimes 
of 10-50 Myr. Dark gaps between the bright rings indicate an absence of dust. Here we show 
that bright, optically thick rings can shadow the disk and produce apparent dark gaps in the disk. 
The rings and gaps produced in our models are comparable in size and surface brightness to those 
observed in some debris disks, such as HR 4796 A. We also derive predicted lifetimes, luminosities, 
and masses for collisional cascades induced by embedded planets or a stellar flyby. Although our 
derived lifetimes are a factor of ~ 2 longer than observed, our results for dust luminosities and 
masses agree with observations (e.g. Greaves & Wyatt 2003; Habing et al. 2001; Wyatt, Dent, & 
Greaves 2003; Decin et al. 2003). 

We outline the model in §2, describe the calculations in §3, and conclude with a brief discussion 
in S4. 



2. THE MODEL 

Kenyon & Bromley (2001, 2002a) describe our multiannulus numerical model for planetesimal 
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growth^ Kenyon &; Luu (1998, 1999), Kenyon (2002), and Kenyon & Bromley (2001, 2002a) 
compare results with analytical and numerical calculations. Briefly, we adopt the Safronov (1969) 
statistical approach to calculate the collisional evolution of an ensemble of planetesimals in orbit 
around a star of mass M (see also Spaute et al. 1991; Weidenschilling et al. 1997). The model grid 
contains N concentric annuli with widths 6ai centered at heliocentric distances Oj. Calculations 
begin with a differential mass distribution n{mik) of bodies with horizontal and vertical velocities 
hik{t) and Vik{t) relative to a circular orbit. The horizontal velocity is related to the orbital 
eccentricity, e^^,(t) = 1.6 {hik{t)/VK,i)'^ , where VK,i is the circular orbital velocity in annulus i. The 
orbital inclination depends on the vertical velocity, i'^i^{t) = sm~^ (2(vik{t) /VK,i)'^) ■ 

The mass and velocity distributions evolve in time due to inelastic collisions, drag forces, and 
long-range gravitational forces. We derive collision rates from kinetic theory and use an energy- 
scaling algorithm to choose among possible collision outcomes. We define So as the tensile strength 
of a planetesimal and Eg as the gravitational binding energy per unit mass. If Ec is the center 
of mass collision energy, the ratio Xc = Ec/{Eg + Sq) sets the collision outcome; Xc ^ 1 yields a 
merger with negligible debris, Xc ~ 1 yields a merger with some debris, and Xc ^ 1 yields only 
debris. The ratio of the collision energy Qf to the crushing energy Qc defines the mass of the 
debris, = Qj/Qc- For most calculations, we adopt a crushing energy Qc = 5 x 10^ erg g^^ and 
a range of So appropriate for icy objects at large distances from the central star, 5o ~ 1 to 10^ 
erg g^^ (Greenberg et al. 1984; Kenyon & Luu 1999). Most fragmentation algorithms adopt a 
minimum velocity for fragmentation, Vf (Davis et al. 1985; Wetherill &; Stewart 1993); we choose 
F/ = 1 cm s-^ (Kenyon & Luu 1999). 

We use two algorithms to compute the velocity evolution from long-range gravitational inter- 
actions. In the high velocity regime, the collision velocity exceeds the mutual Hill velocity, 

VH ~ ^ijaij[{mik + mji)/3M^]^/^ , (1) 

where Qij are aij are the average angular velocity and the average heliocentric distance of annulus 
i and annulus j. Statistical solutions to the Fokker-Planck equation then yield accurate estimates 
for the stirring rates (e.g., Stewart Sz Ida 2000). In the low velocity regime, the Fokker-Planck 
equation underestimates the stirring timescales. For some of our calculations, we adopt the Ida 
Sz Makino (1993) fits to low velocity stirring rates derived from n-body simulations. Kenyon 
& Bromley (2001) describe how we match stirring rates in between the low and high velocity 
regimes. For other calculations, we adopt the Ohtsuki, Stewart, & Ida (2002) stirring rates derived 
from fits to the Stewart & Ida (2000) rates and n-body calculations. In most situations, the 
Ohtsuki, Stewart, & Ida (2002) rates are 20% to 30% smaller than the Kenyon & Bromley (2001) 
rates. The Ohtsuki, Stewart, & Ida (2002) rates also agree better with the analytic formulae of 
Goldreich, Lithwick, & Sari (2002) than the Kenyon &; Bromley (2001) rates. 

Gas drag circularizes the orbits of all mass batches and also removes material from each batch. 



^Thebault, Augereau, & Beust (2003) describe a single annulus model in the context of the /? Pic disk. 
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We adopt a simple nebular model with gas surface density T,g{a,t) = Sgoa~^^^e~*/*9 and scale 
height Hg{a) = HQ^a/aoY'^^^ (Kcnyon & Hartmann 1987) to compute the gas volume density 
Pg. To approximate gas removal on a time tg, the gas density declines exponentially with time. 
We set tg = 10 Myr and adopt the Adachi et al. (1976) formalism to compute inward drift and 
velocity damping from gas drag (see Weidenschilling 1977b, for another discussion of gas drag). In 
calculations at 30-150 AU, particle losses from gas drag are small, ~ 1% or less of the initial mass. 
Velocity damping is negligible at late stages when viscous stirring dominates the velocity evolution 
of small bodies. 

The initial conditions for these calculations are appropriate for a disk with an age of ^ 1-10 
Myr. We consider systems of N annuh in disks with Oj = 30-150 AU and Sai/ui = 0.01-0.025. The 
disk is composed of small planetesimals with radii of ~ 1-1000 m (see below). The particles have 
an initial mass distribution niirruk) in each annulus, with a mass ratio 5 = rriik+i/mik between 
adjacent bins. These objects begin with eccentricity eo and inclination zq = eo/2. Most of our 
models have eo independent of a,; some models have a constant initial horizontal velocity in each 

1/2 

annulus, oc . We assume a power law variation of the initial surface density of solid material 
with heliocentric distance, Sj = Sq (cj/ao)^^/^. Models with Sq ~ 0.1-0.2 g cm~^ at ag = 30 AU 
have a mass in icy solids comparable to the minimum mass solar nebula (Weidenschilling 1977a; 
Hayashi 1981). Observed disk masses for 1-10 Myr old stars range from ~ 10% to ~ 1000% of 
the minimum mass solar nebula (Osterloh &; Beckwith 1995; Natta et al. 1997; Wyatt, Dent, & 
Greaves 2003). 

Our calculations follow the time evolution of the mass and velocity distributions of objects 
with a range of radii, r^k = rmin to rj^ = Vmax- The upper limit Vmax is always larger than the 
largest object in each annulus. To save computer time in our main calculation, we do not consider 

small objects which do not affect significantly the dynamics and growth of larger objects, rmin = 
10-100 cm. Erosive collisions produce objects with rn. < rmin which arc 'lost' to the model grid. 
Lost objects are more likely to be ground down into smaller objects than to collide with larger 
objects in the grid. 

To estimate the amount of dusty debris produced by planet formation, we perform a second 
calculation. Each main calculation yields Mi{t), the amount of mass lost to the grid per annulus per 
timestep. The total amount of mass lost from the planetesimal grid is M = Yli=i -^i(*)- The debris 
has a known size distribution, n'^j = n'^^aj^ , where /3 is a constant. The normalization constant n^Q 
depends only on (3 and M(t), which we derive at each timestep in the main calculation. To evolve 
the dust distribution in time, we use a simple collision algorithm. The optical depth r of the dust 
follows from integrals over the size distribution in each annulus. The optical depth and a radiative 
transfer solution then yield the luminosity and radial surface brightness of the dust as a function 
of time. The appendix describes the collision algorithm and the calculations for the optical depth 
and the dust luminosity. 
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3. CALCULATIONS 

3.1. Planet Formation in a Quiescent Disk 

3.1.1. Wetherill & Stewart Fragmentation 

We begin with calculations in a quiescent disk surrounding a star with a mass of 3 Mq and 
a luminosity of 50 Lq. The disk consists of 64 annuli with Aai/ai = 0.025 and extends from 30 
AU to 150 AU. The initial distribution of icy planetesimals has sizes of 0.5-1000 m with 6 = 2 
and equal mass per mass bin. The initial eccentricity is eo = 10~^ for all objects in all annuli; the 
initial inclination is iq = eo/2. These initial values provide a rough equilibrium between collisional 
damping and viscous stirring at t = 0. Viscous stirring dominates at smaller eo; collisional damping 
dominates at larger eo- The particles have mass density pd = 1.5 g cm^"^. The initial surface density 
in solid objects is comparable to the minimum mass solar nebula, with Sq = 0-18 g cm"^ (ao/30 
AU)-3/2. The total mass in solids, 93.8 M^, is similar to the median dust mass observed in disks 
around young 0.5-3 Mq stars (Wyatt, Dent, & Greaves 2003). We adopt the Adachi et al. (1976) 
gas drag formalism, the Ohtsuki, Stewart, & Ida (2002) stirring rates, and the Wetherill & Stewart 
(1993) fragmentation algorithm with Qc = 5 x 10^ erg g~^ and Sq = 10^ erg g~^. These standard 
choices for Qc and So are appropriate for icy planetesimals with bulk properties similar to terrestrial 
snow (Greenberg et al. 1984). 

Icy planet formation in the outer regions of a quiescent disk has three stages (Kenyon & Luu 
1999; Kenyon 2002). Small planetesimals with ?a 1-1000 m and eo < 10~^ grow slowly. At 
30-37 AU, it takes ~ 0.3 Myr for the largest objects to grow from ~ 1 km to rj ~ 10 km. 
Slow, orderly growth ends when the gravitational cross-sections of the largest objects exceed their 
geometric cross-sections. Runaway growth begins. The largest objects take ~ 2 Myr to reach sizes 
of 100 km and another 15 Myr to reach sizes of 1000 km. As the largest objects grow rapidly, 
dynamical friction and viscous stirring increase the eccentricities of the smallest objects. Collisions 
between these small objects then begin to produce more and more debris. As the largest objects 
grow to sizes of 2000-3000 km, a collisional cascade reduces substantially the mass in the smallest 
objects. During this period of 'oligarchic growth (Kokubo & Ida 1998),' the system reaches a 
rough equilibrium where the largest objects contain most of the remaining mass. 

Figure 1 illustrates the time evolution of the cumulative mimbcr and eccentricity distributions 
at 30-37 AU in the inner disk. The initial distributions are Nc oc with qq = 3 and eq = 10~^ 
for all objects. During slow growth, the velocity distribution develops two features (Figure 1; right 
panel). For the largest objects, dynamical friction dominates viscous stirring and produces a power 
law velocity distribution. At smaller radii, viscous stirring produces a fiat velocity distribution with 
efiat/emin ~ 0.1 (Goldreich, Lithwick, &; Sari 2002). The transition between these two regimes 
occurs at the particle radius where most of the mass is concentrated. During runaway growth, this 
transition moves to larger rf, the shape of the velocity distribution does not change. 

The mass distribution consists of two power laws, with a transition at ~ 1 km (Figure 1; left 
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panel). The largest objects always follow a power law with a; ~ 3. At the start of the slow growth 
phase, collisions produce mergers and little debris. The power law exponent for the small bodies 
thus decreases from the initial ag ~ 3 to ~ 1-2. During runaway growth, debris from collisions 
adds mass to the smallest mass bins; ag slowly increases to the standard collisional exponent, ~ 
2.5 (Dohnanyi 1969; Williams &; Wetherill 1994). Although collisions between small objects reduce 
the total mass contained in small bodies, the slope of the mass distribution is roughly constant, 
with as ~ 2.5. 

Figure 2 shows the evolution of the cumulative surface density in the inner disk. During the first 
1 Myr of evolution, slow growth concentrates most of the mass in 1-10 km objects. Because most 
collisions yield mergers with little debris, the total surface density is roughly constant. Runaway 
growth concentrates more mass in the largest objects. After 10-20 Myr, large objects with radii 
of 100-1000 km contain roughly 10% of the initial mass in the inner disk. Once oligarchic growth 
begins, collisions between small objects produce substantial debris instead of mergers. The surface 
density begins to decline. The largest objects continue to grow and contain an ever-larger fraction 
of the total mass. After 1-2 Gyr, disruptive collisions between small objects reduce the surface 
density by nearly an order of magnitude. Large objects with > 100 km contain ^ 10% of the 
initial mass and ~ 70% of the final mass of solid material at 30-37 AU. 

The timescale for planet growth is a strong function of heliocentric distance. Because the 
collision time is proportional to the surface density and the orbital period, the growth time scales 
with a as t (X P/T, oc (Lissauer 1987; Kenyon &; Luu 1998). Planetesimals thus grow to large 
sizes first in the inner disk. In this model, the ratio of timescales for planets to grow to 1000 km 
at 135 AU and at 33 AU is ~ 80. This result is close to the predicted ratio of ~ # = 64. 

Despite differences in timescales, planet formation proceeds to the same endpoint in all annuli. 
Slow growth, runaway growth, and oligarchic growth produce large objects with ~ 100 km to 
~ 3000 km which contain ~ 5%-10% of the initial mass. These large objects stir up the orbits of 
the leftover planetesimals and initiate a collisional cascade. These disruptive collisions reduce the 
surface density of small objects by roughly an order of magnitude at all a. 

Erosive collisions also lead to structure in the disk (Figure 3). During the slow growth stages, 
most collisions yield mergers. The production rate M of objects with < 1 m is small and 
declines smoothly with increasing heliocentric distance. As planetesimals merge to form planets, 
the collisional cascade produces more and more debris. At the inner edge of the disk, M increases 
by more than 3 orders of magnitude in ~ 50 Myr. The collisional cascade rapidly depletes the inner 
disk of 1-1000 m bodies; M declines. As large planets begin to form in the outer disk, the peak in 
M moves to larger heliocentric radii. This peak moves from Oj ~ 30 AU at t = 80 Myr, to Oj ~ 50 
AU at 400 Myr, and to - 100 AU at t = 2-3 Gyr. 
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3.1.2. Davis et al. Fragmentation 

Calculations with the Davis et al. (1985) fragmentation algorithm allow us to measure the 
sensitivity of the results to the bulk properties of the planetesimals. In the Davis et al. (1985) 
approach, the ejecta from a collision receive a fixed fraction /ke of the center-of-mass collision 
energy. We follow Kenyon & Luu (1999) and adopt Jke = 0.05 and Sq = 1-10^ erg g~^ to 
simulate a wide range of planetesimal bulk properties (see also Greenberg et al. 1984; Davis et al. 
1985). 

Planets form faster with the Davis et al. (1985) fragmentation algorithm (Figure 4; see also 
Kenyon & Luu 1999). During the slow growth phase, collisions that produce some debris tend 
to reduce e and i for the largest planetesimals and increase e and i for the smallest planetesimals. 
Cooling of the large planetesimals increases gravitational focusing factors and speeds the growth of 
the largest objects during runaway growth. As a result, it takes only ^ 15 Myr for planetesimals 
to grow to sizes of 1000 km with the Davis et al. (1985) fragmentation algorithm and Sq = 10^ 
erg g~^, compared to ~ 20 Myr with the Wetherill &; Stewart (1993) algorithm. 

With both fragmentation algorithms, coUisional cascades remove nearly all of the mass from 
a planetesimal disk. Figure 5 compares derived values of M for models with a range of Sq. In all 
models, significant debris production begins when large objects with > 100 km stir small objects 
to the disruption velocity, Vd ^ 50 m (5o/10*^ erg g~^)^/^ (Davis et al. 1985; Kenyon &; Bromley 
2001). Calculations with small 5*0 thus produce more debris at earlier times than calculations with 
large Sq. Planetesimals are easier to fragment, but harder to disrupt, in the Davis et al. (1985) 
fragmentation algorithm than in the Wetherill & Stewart (1993) algorithm. At fixed Sq, models 
with Davis et al. (1985) fragmentation produce more debris at early times, and somewhat less 
debris at later times, than Wetherill & Stewart (1993) models. 

Despite these differences, all calculations have several common features. Shortly after the 
debris production rate reaches a maximum, the decline in M follows a power law with M oc t^^. 
At t ~ 100 Myr, all models converge to M ~ 5 x lO^'' g yr~^. Over a 6 order of magnitude range 
in ^0, the dispersion in M at t 100 Myr is smaller than a factor of two. This common M leads 
to a standard dust mass and disk luminosity at t ~ 10-100 Myr (see below). 

Figure 6 compares final cumulative surface density distributions for the inner eight annuli of 
several models. For So = 10^ erg g~^, 1 Gyr of fragmentation reduces the surface density of solid 
material in the disk by roughly an order of magnitude. For Sq = 1 erg g~^, the cumulative surface 
density is roughly a factor of 30 smaller than the initial surface density. Despite a large range in 
5*0, the range in the final cumulative surface density is small. In these models, large bodies with 
> 100 km contain ~ 1% to 10% of the initial mass and ~ 30% to 80% of the final mass in solid 
material. 

Stochastic processes are an important feature of the collisional cascade. During runaway 
growth, random fluctuations in the collision rate can produce a single large body which grows 
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much more rapidly than bodies in neighboring annuh. By robbing other bodies of material, this 

runaway body slows the growth of large objects nearby. Stirring by the runaway object leads to 
more dust on shorter timescales compared to calculations without a single runaway. Single runaways 
occur in ~ 20% to 25% of test calculations with identical initial conditions. These runaways appear 
to occur more often in models with smaller Sq. 

The general results of the calculations are insensitive to changes in the initial conditions. For 
initial eccentricity cq < 10~^, all calculations yield growth by mergers instead of rapid disruption 
(see Kenyon Sz Bromley 2002a). In models with cq ~ 10~^ to lO"^, coUisional damping reduces e 
until runaway growth begins to produce large objects. Compared to calculations with bq = 10~^, 
growth times are factors of 2-3 longer (see also Kenyon & Luu 1999). Variations in the initial 
surface density, the mass density, or the size distribution of planetesimals change the timescale 
but not the outcome of the evolution. Growth times scale with the initial radius rjo of the largest 
object in the grid. For calculations with Vmin = 0.5-1.0 m, models with r^o ~ 10 m to 1 km have 
indistinguishable growth times when the slope of the initial power law mass distribution is qo < 
4—5. Because the timescale for viscous stirring is shorter than the growth timescale for > 10 
km, models with a substantial fraction of the initial mass in 10 km or larger objects take 2-3 times 
longer to reach runaway growth. Growth rates are more sensitive to the initial mass in solids and 
the mass density of individual planetesimals, t oc T,^^ and t (x pg "^^^ (see Kenyon &: Luu 1999). 

The results are also insensitive to the details of the gas drag and gravitational stirring algo- 
rithms. Because gas drag removes only a few per cent of the initial mass in the disk, changes to 
the damping or drag algorithm do not change the results. The stirring algorithm affects timescales 
for the collisional cascade. Calculations with the Ida & Makino (1993) fits for low velocity stirring 
yield smaller stirring rates than the Ohtsuki et al. (2002) rates. Smaller stirring rates delay the 
collisional cascade but yield similar debris production. For the range of stirring rates we use in our 
calculations, the delay in the collisional cascade is ~ 25%. 

3.1.3. Luminosity and Surface Brightness Evolution 

Our coagulation calculations demonstrate that planet formation in the outer regions of a 
planetesimal disk produces copious amounts of dust. The formation of 1000-3000 km bodies leads 
to a collisional cascade with M ~ 10^^ g yr~^ to lO^'^ g yr~^ at 30-100 AU. For dust survival times 
of ^ 1 Myr at 30-100 AU, these M's imply dust masses of 10~^M^ to lO'^M®. This range is 
comparable to dust masses derived from observations of debris disks (Backman &; Parcscc 1993; 
Lagrange et al. 2000; Zuckerman 2001; Wyatt, Dent, & Greaves 2003). In our simulations, 
collisional cascades yield large M's in narrow ranges of disk radii. The variation in M across the 
face of the disk (Figure 3) suggests a ring of debris that expands in radius as the formation of 
1000-3000 km objects moves to larger disk radii. The apparent dimensions of the ring, Sri ~ 20 
AU at ri ~ 50 AU and Sri ~ 30 AU at ^ 100 AU, are comparable to the dimensions, Sr/r 
0.1-0.2, of the rings in HR 4796A and Vega (Jayawardhana et al. 1998; Koerner et al. 1998; 
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Wilner et al. 2002). 

To quantify comparisons between our models and observations, we calculate the evolution of 
particle numbers for 'dust' with < 1 m. The calculation assumes that all small particles have 
a scale height comparable to the scale height of the 1 m objects in the planetesimal calculation 
and follows the formation and removal of dust grains by collisions and radiative processes. Gas 
drag is not included. The evolution of particle numbers yields the radial optical depth through the 
disk, which we use to derive the time evolution of the radial surface brightness I and the fraction 
of stellar radiation Ld/Lq absorbed and scattered by the disk. This approach complements more 
detailed calculations of dust in a gaseous disk (e.g., Takeuchi k. Artymowicz 2001). Takeuchi & 
Artymowicz (2001) solve the equation of motion for grains in a gas disk irradiated by a central 
star but do not include dust formation by collisions. We derive dust production rates from the 
planetesimal evolution code but do not consider coupling between gas and dust in the disk. The 
appendix describes the details of our calculation. 

Figure 7 shows the evolution of dust mass in several planetesimal disks. We divide the dust 
into small grains with 1 jiia < < 1 mm and large grains with 1 mm < < 1 m. The small 
grains have little mass but produce most of the optical depth; the large grains contain most of 
the mass but have small optical depth. At t = 0, all disks have no mass in small or large grains. 
During the slow growth phase, collisions produce modest amounts of debris. The dust mass grows 
linearly in time. Runaway growth concentrates more and more mass into large objects with > 
100 km, which stir the leftover planetesimals to the disruption velocity. Once runaway growth 
begins in the innermost disk annuli, it takes only 3-10 Myr for the dust mass to grow by 4-6 orders 
of magnitude. The rapid rise in dust mass begins with the formation of objects with > 100 
km. When the largest objects have ~ 1000 km at ~ 30 AU, the dust mass reaches a plateau. 
The dust mass is roughly constant until large objects start to form at the outer edge of the disk. 
Poynting-Robertson drag and radiation pressure then rapidly remove dust from the disk. At late 
times, the decline in dust mass is oc with n ~ 1-2. 

The mass in small or large grains is remarkably independent of the bulk properties of the 
planetesimals or the initial conditions in the planetesimal disk. As planet formation propagates 
through the disk, the dust mass is roughly constant in time. Because planetesimals with small 
Sq are easier to fragment than planetesimals with large ^o, models with Sq = 1 erg g~^ reach 
this plateau more rapidly (~ 1 Myr) than models with = 10^ erg g~^ (30-50 Myr). Despite a 
six order of magnitude range in So, the range in dust mass is a factor of two for similar starting 
conditions. Minimum mass solar nebula models with an initial mass of ~ 100 in solids at 
30-150 AU yield a mass in small grains of Mg ~ 0.01-0.02 M®; the mass in large grains is Mi ~ 
0.1-0.25 M®. The dust mass is equally insensitive to large ranges in cq or pd- The dust mass is 
roughly proportional to the initial mass of solids in the disk, Mj, Mi oc Mq. 

The dust masses derived from the planetesimal model are comparable to those observed in 
debris disk systems. The masses in small grains are close to the minimum mass needed to produce 
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an observable mid-IR or far-IR excess of radiation for a nearby main sequence star (e.g., Wood 
et al. 2002). We derive Ms ^ 0.01-0.02 M© in the plateau and Ms > lO'^M® at late stages 
of the evolution. Backman & Paresce (1993) quote minimum masses of ~ 0.001-0.01 for a 
Lyr, a PsA, and /3 Pic. For a sample of remnant disks around nearby main sequence stars, Habing 
et al. (2001) quote minimum dust masses of 10-^-10"^ (see also Greaves Sz Wyatt 2003; 
Wyatt, Dent, & Greaves 2003; Decin et al. 2003). 

During the coUisional cascade, large dust masses result in luminous disks. Figure 8 illustrates 
the time evolution of the relative disk luminosity Ld/Lq for models with different Sq. The left panel 
shows luminosity evolution for models without radiation pressure on smaU grains. The right panel 
shows the luminosity evolution when radiation pressure removes small grains with < 1 /Ltm on 
the local dynamical timescale. Independent of 5*0, models with radiation pressure yield maximum 
disk luminosities of Ld/Lq ~ 10^^. Because radiation pressure is unimportant for large Sq, the 
maximum disk luminosity is independent of radiation forces for > 10^ erg g^^. In models with 
small So, very small grains ejected by radiation pressure contribute most of the optical depth in 
the disk. Thus, models without radiation pressure yield lower disk luminosities for Sq < 10^ erg 



Although the magnitude of the disk luminosity is sensitive to the treatment of radiation pres- 
sure, the form of the luminosity evolution is independent of Sq and other parameters in the cal- 
culations. All models have a relatively rapid rise in Lu/Lq followed by a longer decline. While 
the mass in dust is relatively constant in time, the luminosity follows a shallow power law decline, 
Ld/Lq oc t~™', with m 1. As planet formation propagates out through the planetesimal disk, dust 
forms at larger and larger distances from the central star. Because dust in the outer disk intercepts 
less radiation from the central star than dust in the inner disk, the dust luminosity declines with 
time. Once Poynting-Robertson drag removes material from the disk, the luminosity declines more 
rapidly with a power law index closer to m = 2. 

CoUisional cascades also yield a standard maximum luminosity for the outer part of the disk. 
The luminosity is insensitive to the bulk properties of planetesimals, fxE, Pd, Qc, Sq, and Vf, and 
many of the initial conditions, cq, qq, and ViQ. The luminosity is sensitive to the initial mass in 
planetesimals Mq at 30-150 AU: 



This luminosity is comparable to the maximum luminosity of known debris disks, such as (3 Pic and 
HR 4796A. Because the timescale to reach the maximum luminosity depends on the initial mass in 
planetesimals, we can write the dust luminosity as a function of time. 



g 



.-1 




(2) 




(3) 



where to ^ 30 Myr (Mo/100 M^)''^. The coefficient for to ranges from 1-3 Myr for Sq = 1 erg 
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g ^ to ~ 100 Myr for = 10^ erg g ^. Dominik &; Decin (2003) derive similar relations from an 
analytical model for the collisional cascade. 

Finally, all planetesimal calculations produce axis3niimetric structures within the disk. Figure 9 
shows two sets of relative surface brightness distributions for models with different Sq. For = 10^ 

7/2 

erg g~ (top panel), the initial surface brightness profile in curve (a) is a power law, I/Iq cc a- 
During the slow growth phase, the surface area per unit mass of the planetesimals drops and the 
surface brightness fades. Because planets grow faster in the inner disk, the surface brightness of 
the inner disk fades more rapidly than the brightness in the outer disk. The surface brightness 
reaches a minimum at t = 4 Myr; I/Iq oc a~'^^^- Once runaway growth begins, it takes only 7 Myr 
for the surface brightness at the inner edge of the disk to reach a maximum roughly two orders of 
magnitude brighter than the initial planetesimal disk. As planet formation propagates through the 
disk, a bright ring of dust emission moves outward. This ring highlights the region of maximum 
dust production and signals the presence of at least one planet with a radius of 1000 km or larger. 
It takes ~ 60 Myr for this ring to reach a, ~ 50 AU and another 330 Myr to reach 80 AU. If we 
define the width Sa/a of the ring by the radius where the surface brightness drops to half of the 
maximum, Sai/ai ~ 0.2 at = 40 AU and Sai/ai ~ 0.15 at = 80 AU. At t = 400 Myr. the 
inner disk is nearly two orders of magnitude fainter than at t = 0. After 2 Gyr, planet formation 
reaches the outer edge of the disk and the surface brightness of the entire disk fades. As the disk 
fades, the surface brightness rises with heliocentric distance, with I/Iq oc and p 0-2. 

CalcTilations with smaller Sq produce shadowed disks (Figure 9; bottom panel). During the 
early stages of models with Sq < 10^ erg g~^, the disk surface brightness is a power law and fades 
slowly with time (curves (a) and (b)). Once runaway growth begins, the inner disk brightens by 
2-3 orders of magnitude. The dust in the inner disk is optically thick, r 3-10 at Oj pa 30 AU, 
and shadows material in the outer disk. Shadowing produces a pronounced minimum in the surface 
brightness, which propagates outwards as planets form at progressively larger disk radii. Although 
this shadow resembles the dark gaps cleared of material by a large planet, it is not an absence of 
dust. The shadow is a region where light from the central star does not penetrate and is therefore 
of much lower surface brightness than surrounding bright regions with comparable amounts of dust. 
In models with 5*0 = 10^ erg g^^, it takes only ~ 2 Myr to form a bright inner ring and a dark 
shadow at ~ 30-40 AU. This structure moves to ~ 50-60 AU at 20 Myr and to ~ 80-100 
AU at 100 Myr. The radial extent of the dark gap is 5a/ai ~ 0.1. After 500 Myr to 1 Gyr, planet 
formation in the outer disk starts to produce dust and the entire disk begins to fade. 

Figures 10-11 show a montage of color snapshots of the disk surface brightness. In Figure 10, 
models with 5o = 10^ erg g~^ produce a bright ring with a dark shadow outside it. As the ring and 
shadow propagate out through the disk, they fade relative to the surface brightness of surrounding 
disk material. In Figure 11, models with = 10^ erg g~^ produce a bright ring that propagates 
out through the disk. The ring is brightest where planets with radii of ~ 1000-3000 km form. 
Inside the ring, planet formation has saturated with the formation of 3000 km objects. Outside 
the ring, large planets have yet to form. 
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In the animations of Figures 10-11 included in the electronic version of this paper, planet 
formation appears as a succession of waves flowing outward in the disk (see also Kenyon & Bromley 
2002b). Slow growth from 1 km to ~ 100 km produces little dust and concentrates more and more 
material into objects with a smaller geometric cross-section per unit mass. Thus, a dark wave 
moves out through the disk. Dust formed during runaway growth lies in bright rings which appear 
as a bright wave moving out through the disk. Once planet formation is complete, a dark wave - 
which heralds the disappearance of dust - propagates out through the disk. During this phase, the 
largest bodies in the disk may coalesce to form larger planets. 

3.2. Planet Formation After a Stellar Flyby 

Planetesimal calculations demonstrate that embedded planets with > 1000 km can produce 
bright rings and dark gaps in a quiet planetesimal disk. Stellar flybys can also produce rings, gaps, 
and possibly other structures in a planetesimal disk (Larwood 1997; Mouillet et al. 1997; Kalas 
et al. 2000; Ida et al. 2000; Kobayashi & Ida 2001; Larwood & Kalas 2001). The lifetimes for 
rings and gaps produced by a flyby are much shorter, < 1 Myr, than the lifetimes of rings and gaps 
produced by planets, ~ 10-100 Myr. Because field stars arc less likely to interact with passing stars 
than stars in young associations, flybys are more relevant to understanding the structures observed 
in young debris disks than in old debris disks (e.g. Ida ct al. 2000; Kalas et al. 2000; Kenyon &: 
Bromley 2002a). Here we consider whether planet formation can regenerate bright rings and gaps 
following a stellar flyby. 

Because our code does not currently follow the trajectories of individual objects, we assume 
that the close passage of a low mass star at a > 600 AU instantaneously raises the eccentricities of 
planetesimals. This approximation is valid if the encounter between the passing star and the disk is 

short compared to the collisional damping time. Passing stars not bound to the central star of the 
disk satisfy this requirement (Larwood 1997; Mouillet et al. 1997). We adopt a functional form 
for the imposed eccentricity that satisfies constraints derived from ri-body calculations (Larwood 
1997; Mouillet et al. 1997; Kobayashi & Ida 2001). 

Kenyon & Bromley (2002a) describe the early evolution of a planetesimal disk following a 
moderately close stellar flyby. Large perturbations with eo > 0.05 lead to complete disruption of 
nearly all planetesimals (Kenyon & Bromley 2002a, and references therein). When a flyby produces 
a modest perturbation with eo ^ 0.03-0.04 at 30-150 AU, two body collisions produce substantial 
amounts of dust and damp planetesimal velocities on short timcscalcs. Once planetesimals have e < 
0.01, collisions produce growth instead of disruption. Continued damping and growth eventually 
produces planets. 

To illustrate planet formation in a planetesimal disk following a stellar flyby, we consider in 
detail a disk model with eg = 0.02 (cj/ao)^''^, Sq = 0.02 g cm^^ at gq = 30 AU, and a total mass 
in solids of 9.3 M^. This disk mass is at the lower end of the observed range of dust masses for 
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0.5-3 stars (Wyatt, Dent, & Greaves 2003). The disk consists of 64 annuli with Aai/ai = 
0.025 and extends from 30 AU to 150 AU. The initial distribution of planetesimals has sizes of 10 
cm to 10 m with 6 = 2 and equal mass per mass bin. The particles have mass density pd = 1.5 g 
cm~^. We adopt the Kenyon & Bromley (2001) stirring algorithm and the Wetherill &: Stewart 
(1993) fragmentation algorithm with Qf. = 5 x 10^ erg g"-*^ and Sq = 2 x 10^ erg g~^. Kenyon & 
Bromley (2002a) describe several aspects of the early evolution of this model. 

Planet growth following a stellar flyby has four stages. After the flyby, planetesimals damp 
quickly and begin to form larger objects. This evolution starts at the inner edge of the disk and 
slowly propagates outward. At all heliocentric distances, the 50-200 cm planetesimals damp first. 
This damping produce a V-shaped eccentricity distribution in < 1 Myr at 30-37 AU (Figure 12). 
After 2-4 Myr, larger planetesimals damp and begin to grow slowly. As objects grow to sizes of 
0.1-1 km, dynamical friction dominates viscous stirring and reduces the eccentricities of the largest 
objects. 

When particles reach sizes of 1 km or larger, the growth and stirring rates increase rapidly 
(Figure 13). Damping of the eccentricities of the largest objects leads to runaway growth. At 30-37 
AU, it takes 40 Myr for the largest objects to reach sizes of 10 km and another 20 Myr for large 
objects to reach sizes of 100 km. These timescales are a factor of ~ 2 larger than growth times 
for icy planets in a quiet disk. Viscous stirring by the largest bodies increases the eccentricities 
of the smallest objects. The velocity evolution quickly reaches an equilibrium between the largest 
bodies with radii of 100 km or larger and 1-10 km bodies containing most of the mass. This ratio 
is roughly eiarge/^ small ~ 0.03-0.1 (Goldreich, Lithwick, &; Sari 2002). As the largest bodies grow 
to sizes of 1000 km, they stir up the smallest bodies to this equilibrium eccentricity ratio. 

All processes take longer farther out in the disk. The timescales for collisions, damping, and 
stirring grow with heliocentric distance. The collision time is t oc P/S oc a^, which implies a 
damping timescale of ~ 0.3 Myr at 30-37 AU and ~ 25 Myr at 130-150 AU. Because our model 

1/2 . 

has eo oc a/ , the outer disk damps more slowly and loses more material than the inner disk. This 
extra mass loss slows the rate of planet formation by another 20% to 30% relative to the inner disk. 
Thus, planet formation takes ~ 200 times longer at 150 AU than at 30 AU. 

Throughout the evolution, the cumulative size distribution has three main components (Figure 
14). As collisions damp particle eccentricities, fragmentation dominates growth by mergers. The 
size distribution for the small objects follows a power law with w 2.5 for r < 0.1-1 km. Once 
mergers become important, the size distribution develops a merger component, with a; pa 3 for r > 
1-10 km. At intermediate sizes, 0.1 km < r < 10 km, the size distribution has a pronounced hump 
containing most of the mass. As planets grow in the disk, the position of this hump moves to larger 
radii. 

The right panel of Figure 13 illustrates some features of the collisional cascade induced by 
planet formation. During the early evolution of this model, erosive collisions remove material from 
the disk. The cumulative surface density declines by ~ 30% at 30-37 AU and 40% to 50% at 125 — 
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150 AU. As mergers produce 1-100 km objects, the cumulative surface density is roughly constant 
in time. Once 100-1000 km objects begin to form, viscous stirring increases the eccentricities of 
the leftover planctcsimals (Figure 12). Continued stirring leads to the coUisional cascade, where 
0.1-1 km planetesimals are slowly ground into smaller and smaller objects. Because 1-10 km 
planetesimals are too strong, the hump in the size distribution shifts from ~ 1 km to ~ 10 km. 
After 2 Gyr, the collisional cascade produces a factor of two decline in the cumulative surface 
density. 

Planet growth and the collisional cascade depend on the initial mass in planetesimals. More 
massive disks have shorter collision times, with tc oc M^^ (e.g., Lissauer 1987; Kenyon Sz Luu 
1998, 1999). In a disk with a mass in solids comparable to the minimum mass solar nebula, planets 
grow a factor of ten faster than in the model in Figures 11-13. Figure 14 shows that the character 
of the evolution is not sensitive to the initial mass. The cumulative number distribution rapidly 
develops a smooth power law for r < 0.1 km, with ag ~ 2.5. For r > 10 km, mergers lead to a 
steeper power law with a; f« 3 on longer timescales. After ~ 30 Myr, the largest objects have r ^ 
1000 km at 30-37 AU. The collisional cascade begins and reduces the surface density by a factor 
of 3 over the next 300 Myr. After 3 Gyr, erosive collisions reduce the surface density at the inner 
edge of the disk by another factor of three. The largest objects then have radii of ^ 2000-3000 km; 
most of the mass is in objects with r ~ 10-100 km. 

To test the robustness of these results, we calculated models with a variety of initial conditions. 
We changed 5*0, eo, and the initial variation of e with a. Besides the initial mass in planetesimals, 
the ratio eo/So is the important parameter in these calculations. Gompared to our baseline models 
with eo/So = 10~^, models with eo/So > 10~^ lose a larger fraction of their initial mass and take 
longer to form planets. When eo/So>3 x 10~*, erosive collisions remove almost all planetesimals 
before collisional damping becomes effective. These calculations do not form planets on reasonable 
timescales, < 1 Gyr. For calculations with fixed cq/Sq, the mass of the largest planet scales with 
5*0. In our models, collisional cascades begin sooner when Sq is smaller. Because the collisional 
cascade robs material from the largest bodies, smaller 5*0 prevents large objects from growing. 

Calculations with the Davis et al. (1985) fragmentation algorithm yield similar results. Our 
tests show that damping is more efficient in calculations where fx e is large and less efficient when 
/ke is small. In general, icy planetesimals with fxE 0.03-0.10 are harder to disrupt and easier 
to fragment in the Davis et al. (1985) algorithm. Thus, these models tend to produce less debris 
and to form planets on shorter timescales than our baseline models. This difference is small, ^ 
25% to 35%. 

Stellar flybys and planet formation produce copious amounts of dust. Figure 16 compares the 
evolution of the debris production rate M of flyby models with the M evolution for a quiet disk 
model. During the flyby, collisions produce substantial debris on a short timescale. The derived M 
is 6-7 orders of magnitude larger than the initial M for a quiescent disk and 4-5 orders of magnitude 
larger than the largest M for a planet-forming disk. GoUisions damp planetesimal velocities; M 
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declines. As collisions start to favor mergers over debris production, the decline in M accelerates. 
When 10 km objects begin to form in the inner disk, the global M reaches a minimum. The growth 
of 100 km objects leads to a linear increase in M; the growth in M accelerates once 1000 km objects 
begin to form and stir up the leftover planetesimals. Debris production saturates during oligarchic 
growth and then declines as planet formation propagates out through the disk. 

Compared to quiet disk models, planet formation following a stellar flyby produces less dust 
at later times (Figure 16). After a stellar flyby perturbs a planetesimal disk, planets grow more 
slowly; once planet formation begins, the reservoir of leftover planetesimals is smaller. For a modest 
perturbation with cq = 0.02, dust production following a flyby begins 3-10 times later and peaks 
at a factor of two smaller M than dust production in a quiet disk. Because A-type stars evolve 
into red giants on 1 Gyr timescales, stronger perturbations probably prevent planet formation at 
30-150 AU around A-type stars. 

We conclude this section with a brief discussion of Figure 17, which illustrates the evolution 
of the dust luminosity following a stellar flyby. The initial perturbation produces a substantial 
dust luminosity that is 2-4 times larger than the maximum luminosity of a planet-forming disk 
(compare with Figure 8). This large luminosity is short-lived and declines rapidly as collisions cool 
the planetesimal swarm. As objects grow to 10 km sizes in the inner disk, the dust luminosity 
reaches a minimum roughly 3 times more luminous than the luminosity minimum for a quiescent 
planet-forming disk. Dust production associated with planet formation leads to a rise in luminosity 
followed by a linear decline as in the quiet disk models. The peak luminosity in flyby models is 
significantly smaller than in quiet disk models. For similar initial masses, flyby models with eo = 
0.02 are a factor of 3-5 fainter than quiet disk models. The flyby also causes a substantial delay 
in the rise in dust luminosity. Flyby models with eo = 0.02 reach maximum dust luminosity 5-10 
times later in time than quiet disk models. 



3.3. Limitations of the Models 

In previous papers, we have described limitations to multiannulus (Kenyon & Bromley 2001, 
2002a) and single annulus (Kenyon & Luu 1998, 1999) coagulation calculations. Here, we summa- 
rize the most important of these limitations and consider uncertainties in the dust calculation. 

As long as the statistical assumptions underlying the formalism are met, coagulation calcu- 
lations provide a reasonable representation of real collision evolution (Wetherill 1980; Greenberg 
et al. 1984; Davis et al. 1985; Barge & Pellat 1991; Spaute et al. 1991; Lissauer & Stewart 
1993; Wetherill & Stewart 1993; Stern & Colwell 1997; Weidenschilling et al. 1997; Kenyon & 
Luu 1998; Liaba et al. 2001). In our calculations, the spacing of mass bins in an annulus and the 
spacing of annuli in the disk limit the accuracy of the results. Our choice of mass spacing, 5 = 2, 
lengthens the evolution time by 10% to 20% (see Kenyon & Luu 1998, and references therein). 
The radial resolution in our grid of annuli, Aaj/aj = 0.025 also lengthens the evolution time. Tests 
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with different grids suggest a lag of 15% to 25% relative to calculations with very fine grids. Tfius, 
our evolution timescales are ~ 30% to 40% longer than the actual evolution times. 

If N is the number of annuli and M is the number of mass batches, the computation times 
increase roughly as N'^M'^. Because a typical calculation requires 2 4 weeks of cpu time, increasing 
the resolution of a calculation to achieve a more accurate evolution time is prohibitively expensive. 
Improving the accuracy of our evolution times requires practical increases in computing speed, 
which will be achieved with the next generation of high-speed parallel computers. 

The coagulation algorithm begins to break down when binary interactions between large objects 
become important. We reach this limit when ~ 1000-2000 km. At this point, the coagulation 
algorithm underestimates collision and stirring times, and overestimates the evolution time for 
dust clearing and the growth of the largest objects. Tests with a hybrid n-body-coagulation code 
(Bromley & Kenyon 2003) suggest that our pure coagulation results overestimate the dust clearing 
time by less than a factor of two. Although the hybrid code forms larger objects than the pure 
coagulation code, differences in the final mass distribution are small. We plan to report on these 
calculations in a future paper. 

In our implementation, the inherent limitations of the coagulation algorithm have clear ob- 
servational consequences. We overestimate the timescale to produce large planets by ~ 40% and 
the timescale for dust clearing by a factor of ~ 2. Because the collisional cascade must remove 
the same amount of material on a shorter timescale, our estimates for dust mass and luminosity 
are probably low by a factor of ~ 1.5-2.0 compared to ideal calculations with perfect resolution in 
mass and radius. Our bright dust rings are probably 50% larger in radius than those in an ideal 
calculation; our dark gaps are probably smaller by a similar factor. If the actual clearing time for 
dust is shorter, then our radial surface brightness profiles are too shallow during the late stages 
of the evolution. All of the uncertainties in timescales are small compared to the 1-2 order of 
magnitude range in evolution timescales set by the range in initial disk mass (e.g. Wyatt, Dent, & 
Greaves 2003) and the bulk properties of planetesimals. Thus, the observational uncertainties are 
set more by unknown physics than by limits in the calculations. 

Our estimates for the evolution of dust mass, luminosity, and surface brightness rely on several 
assumptions, (i) large bodies do not accrete small grains, (ii) small grains have the same scale height 

—5 /2 

as 1 m objects, and (iii) the size distribution of dust grains is fixed at A'^c oc r • . Large bodies 
probably accrete ~ 1-3% of the total mass converted to dust, which has a negligible impact on the 
evolution of small or large bodies. Stirring by large planets sets the scale height of the smallest 
objects in the planetesimal grid. The amount of stirring does not depend on the mass of the small 
objects. Because gas drag is unimportant at late stages in our calculations, the scale height of 
small objects is fairly independent of their mass (Figures 1, 4, and 8). We estimate a factor of 
^ 2 uncertainty in the scale height of the small objects, which yields factor of two uncertainties 
in the dust luminosity. Dohnanyi (1969) and Williams & Wetherill (1994) show that collisional 
cascades produce power law size distributions with ag ~ 2.5 (see also Tanaka et al. 1996). The 
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uncertainty in this exponent is small, ±0.05 or less for a large range in input parameters. In our 

calculations, significant deviations from this power law can occur at a boundary size, rj,, where 
collisions produce complete disruption for rj < and modest fragmentation for > rj. For the 
bulk properties assumed in our calculations, r;, > 1 m, the minimum size of planetesimals in our 
main calculation. Thus we expect Ug ~ 2.5 for small particles to a high degree of accuracy. This 
result yields a small uncertainty, < 5%, in the derived dust masses and luminosities. 



4. DISCUSSION AND SUMMARY 

Our calculations demonstrate that icy planet formation is the inevitable outcome of coagulation 
in the outer regions of a planetesimal disk (see also Lissauer 1987; Kenyon 2002). In a quiet disk, 
icy planets with r-j !^ 1000-3000 km form on a timescale 

where Timmsn is the surface density of the minimum mass solar nebula. Planets form more slowly 
in a disk perturbed by a stellar flyby. Large perturbations prevent planet formation. In modest 
perturbations with cq ~ 0.02-0.04, planet formation timescales are a factor of 2-4 longer than 
timescales in a quiet disk. 

Icy planet formation produces copious amounts of dust. As planets grow to sizes of 100-1000 
km, gravitational stirring increases the orbital eccentricities and inclinations of leftover planetes- 
imals. High velocity collisions between leftover planetesimals lead to a coUisional cascade, where 
1 km objects are slowly ground to dust. The collisional cascade reduces the surface density in the 
disk by 90%-95%. The reduction in surface density depends weakly on the tensile strength 5*0 of 
1 km planetesimals. Collisions between 'weak' planetesimals with ^ 10^ erg g~^ produce more 
dust than collisions between 'strong' planetesimals with ^ 10^ erg g^^. 

In all of our calculations, dust is confined to bright rings separated by dark gaps. Bright rings 
form along the orbits of large planets, which stir leftover planetesimals up to the disruption velocity. 
Dark gaps are regions with little dust, where planets have not yet formed, or regions shadowed by 
bright, optically thick dusty rings in the inner disk. The bright rings and dark gaps have sizes, 
0.1-0.2 a, comparable to the rings and gaps observed in HR 4796A and other debris disks. 

Dusty rings produced by icy planet formation are observable. The maximum luminosity of 
a dusty disk, equation (2), is comparable to the dust luminosities of known debris disk systems 
(Backman & Paresce 1993; Kalas 1998; Lagrange et al. 2000; Habing et al. 2001; Spangler et 
al. 2001; Greaves & Wyatt 2003; Decin et al. 2003). The decline in dust luminosity follows a 
power law, oc with m k. 1-2 (Dominik & Decin 2003, see also equation (3)). If the initial 
mass in planetesimals is comparable to the mass of a minimum mass solar nebula, the timescale 
for the decline in dust luminosity is 500 Myr to 2 Gyr. More massive disks evolve on shorter 
timescales. The maximum dust luminosity agrees with observed luminosities of debris disks. Our 
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derived lifetimes are longer than the observed lifetimes (see also Greaves & Wyatt 2003; Decin et 
al. 2003). Limitations in our calculations probably lead us to overestimate the decline timescales 
by a factor of ~ 2-3. Given the uncertainties in the initial conditions, this 'correction' brings our 
model lifetimes into reasonable agreement with observations. 

Aside from the initial disk mass and the mass density of planctcsimals, our results are in- 
sensitive to most input parameters. The dust mass - and thus the dust luminosity and surface 
brightness - does not depend on the bulk properties or the initial orbital or mass distributions of 
the planetesimals. The evolution timescales are longer if planetesimals are more dense, stronger, 
and have larger initial orbital eccentricities. Porous, weaker planetesimals in more circular orbits 
grow more rapidly. For fixed initial conditions, the uncertainties in the results are small compared 
to the 1-2 order of magnitude range in the initial disk mass (Wyatt, Dent, & Greaves 2003). 

Together with other studies, our results demonstrate that several physical processes can pro- 
duce rings, gaps, and other structures in a debris disk. If the gas density in the disk is modest, 
coupling between the gas and the dust can produce rings and gaps (Takeuchi & Artymowicz 2001). 
Resonant interactions between dust and distant large planets (Wilner et al. 2002; Kuchner k. Hol- 
man 2003) or small moons (Wyatt et al. 1999) can produce rings or ring arcs in the disk. Large 
planets can clear large gaps in the disk (Lin &; Papaloizou 1979; Goldreich & Tremaine 1980). 
Our models produce bright rings from collisional cascades induced by small planets within the ring. 
Dark gaps arise between small planets in nearby bright rings or by shadowing of the disk by an 
optically thick bright ring. 

In all models, the visibility of dusty structures in a debris disk is related to the collision rate. 
The collision lifetime of a single particle of radius r in an annulus of width Aa at distance a from 
the central star is tc ~ {dn/dt)~^: 

a Aa , , 

tc ~ 2 ' ^ ' 

where $7 is the orbital frequency and n is the number of particles. The luminosity of the annulus 
is roughly Ld/Lq « TH/a « {r/af'n, which yields 

t.~10Myrf^) f-^) (6) 

\ a J VlOOO yry ^ IQ-s J ^ ' 

Because the optical depth and collisional cross-section depend on r^, the collision time is indepen- 
dent of particle size. Most debris disks have relative luminosities Ld/Lq > 10^^, orbital periods 
P < 1000 yr, and relative sizes 6a/a ~ 0.1-0.2. The collision time of 10 Myr or less is shorter 
than the typical stellar age, t > 10-100 Myr. Particle collisions thus play an important role in 
the evolution of all known debris disks. Collisions are most important in the most luminous debris 
disks such as /? Pic and HR 4796A. 

As samples of debris disk systems grow, observations might measure the importance of col- 
lisions. In bright disks, very small grains ejected by radiation pressure produce a large optical 
depth. Systematic variations in the radial scale height with disk luminosity or stellar age might 
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imply some emission from these small grains. Measured disk luminosities inside bright rings might 
constrain the amount of material dragged inwards by Poynting-Robertson drag and thus the mass 
in small grains inside the bright ring. As our calculations become more sophisticated, we plan to 
make detailed predictions for comparisons with observations. 

Our results suggest caution when interpreting bright and dark structures in debris disks (see 
also Takeuchi &: Artymowicz 2001). Dark gaps can indicate a large planet (Lin &: Papaloizou 
1979; Goldreich k, Tremaine 1980), with a mass nip k, 3(Aa/a)^M*. For the measured sizes of 
dark gaps in known debris disks, rUp/M^ ^ 10"^. In our models, shadowing or several bright 
rings can produce dark gaps with Aa/a !^ 0.01-0.1. These structures imply smaller planets with 
nip/Mi, ^ 10~^. In systems with residual gas, interactions between the gas and small dust grains 
can produce a single bright ring (Takeuchi &; Artymowicz 2001) without any planet. 

Observations can plausibly distinguish between these possibilities. Multiple rings, apparently 
observed in (3 Pic (Wahhaj et al. 2003), favor planet formation over dust-gas interactions. Measured 
scale heights or velocity dispersions might distinguish between large or small planets associated with 
a bright ring or a dark gap. Large scale heights with H/a 0.1 are more consistent with Jovian 
mass planets; small scale heights with H/a 10~^ imply smaller planets with masses ranging from 
Pluto to the Earth. 

Finally, multiannulus accretion codes are an important step towards understanding the for- 
mation and evolution of planetary systems (see also Spaute et al. 1991; Weidenschilling et al. 
1997). Our calculations demonstrate that 'complete' disk models with 64 or more annuli provide 
interesting conclusions regarding icy planet formation in a planetesimal disk. As improvements in 
computers allow larger simulations, we plan to extend these calculations into the terrestrial and 
giant planet regimes to improve constraints on the planetesimal theory. These calculations also 
make interesting predictions regarding outcomes of planet formation. Observations with satellites 
- e.g.. SIRTF and JWST - and ground-based telescopes - e.g., OWL and GSMT - will provide 
tests of these predictions and lead to improvements in our understanding of planet formation. 

We acknowledge a generous allotment, ~ 1500 cpu days, of computer time on the Silicon 
Graphics Origin-2000 computers 'Alhcna', 'Castor', and 'Pollux' at the Jet Propulsion laboratory 
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We thank P. Goldreich for clarifying our understanding of gravitational stirring in the low velocity 
limit and for suggesting plots of cumulative surface density. Portions of this project were supported 
by the NASA Astrophysics Theory Program, through grant NAG5-13278. 
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A. APPENDIX 

Our planetesimal calculations derive the evolution of the numbers and orbital elements for 
objects with radii of 1 m and larger as a function of time (Kcnyon & Bromley 2001, 2002a). Detailed 
calculations for the evolution of the smaller particles is prohibitively expensive in computer time. 
To estimate the number of smaller particles, we perform a second, approximate calculation that 
uses the production rate M of small particles derived in the planetesimal calculation^. Throughout 
the later stages of the planetesimal calculations, /ijfc ~ constant for 1 m to 1 km objects. We thus 
assume that the scale height of small particles is identical to the scale height of 1 m objects in the 
planetesimal grid. During a collisional cascade, Dohnanyi (1969) and Williams & Wetherill (1994) 
show that the cumulative number distribution of small objects is a power law, Nq = iVor~", with 
a ~ 2.5. At late times, the 1 m to 1 km objects in our planetesimal calculations reproduce this 
reult. We therefore assume that the dusty debris produced in the planetesimal calculations has 
the same power law. The dust production rate in each annulus Mj thus yields the normalization 
constant Nq of each annulus for each timestep. 

We consider two populations of dust, very small grains and larger grains. On a dynamical time 
scale set by the local circular velocity, vk = (GM/a)^^"^ , radiation pressure ejects very small grains 
with radii between ri and r2- If M^g is the production rate of very small grains in annulus k, the 
space density of very small grains ejected from annulus k is 

Pk = -A — 2 ^-r ' (^1) 

ATrapKk sm Ok 

where 9k is the opening angle defined by the vertical scale height of the very small grain population. 
The volume density pk depends on the size distribution of very small grains. 



47rp, 

Pk = 



Tw^-", (A2) 



where pg is the mass density of an individual grain. We adopt a = 2.5 and assume that grains have 
the same mass density as particles in the planetesimal grid. Solving the two equations for pk yields 
the set of normalization constants raofc- The normalization constants yield the volume density of 
the outflowing wind of very small grains as the sum of material ejected from all interior annuli 



Y.p'^ (A3) 



k=l 



To derive the optical depth of the very small grains, we adopt the geometric optics limit 

Tsi = 2vr / n{r, a)a^dadr, (A4) 



^Thebault, Augcrcau, & Beust (2003) perform a similar calculation in a single annulus and apply their results to 
the inner disk of /3 Pic. 
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where n = no^r " is the space number density. We solve the integral over the size distribution 
exactly and sum the optical depth through the planetesimal grid: 



_ 3(^r2/ri-l) ^ 
87rpgr2(l - \/ri/r2) -^^ 

where 05 is the inner boundary of an annulus centered at o^. 




(A5) 



Although the exact expression for the optical depth is a complicated function of grid variables, 
we derive a simple expression based on characteristic quantities averaged over the grid. For a 
reasonable grain population, ri = 0.01 //m and r2 = 1 /xm. The optical depth is 



M \ /30AU\ /lOkms-^\ /IQ-^X /l.5gcm-3\ AO 



-4 



cm 



When the production rate of very small grains exceeds 10^^ g yr~^, the optical depth approaches 
unity. This rate corresponds to the ejection of 1-2 Earth masses every million years. 

To derive the time evolution of the number density and optical depth for the larger grains, 
we consider three processes. Erosive collisions between particles in the planetesimal grid yield 
large grains with a size distribution Nc = NQr~'^'^. Erosive collisions between these particles 
remove material from the large grain population but maintain the same power law size distribution. 
Poynting-Robertson drag preferentially removes small grains and changes the size distribution. In 
each annulus, we begin with a set of discrete mass bins with minimum radius r2 and maximum 
radius Vmax', fmax is the radius of the smallest particle in the planetesimal grid. In most of our 
calculations, r„j„.^. = 1 m. At t = 0, all of the large grain bins are empty. For each timestep, we 
add particles to the bins using Nc derived from the dust production rates Mj. We compute the 
collision rates and outcomes for the largest mass bin in each annulus and scale this rate to derive 
collision rates and outcomes for smaller dust grains. We add or remove mass from each bin based 
on these rates. Finally, the change in particle number due to Poynting-Robertson drag is 

= \jkinik (A7) 

where Xij^i depends on the size, density, and radiative properties of the grains, the stellar flux, and 
the relative numbers of grains in adjacent annuli (e.g. Burns, Lamy, & Soter 1979). Because these 
properties change slowly with time, we assume Xij^i is constant during each timestep and integrate 
driik/dt exactly to derive the amount of material dragged through the grid. This algorithm yields 
the number of grains in each mass bin in each annulus as a function of time. To derive the optical 
depth of the dust in each annulus, we assume the geometric optics limit and sum the optical depth 
of each mass bin. 

Our simple dust collision algorithm yields the optical depths in planetesimals and planets, 
large dust grains, and very small dust grains in a discrete set of concentric annuli surrounding a 
star. The optical depth of grains dragged out of the innermost annulus by Poynting-Robertson 
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drag is small compared to the optical depth in all of the annuli; we ignore this contribution to the 
optical depth. Tests of the algorithm indicate that the optical depths are accurate to a factor of 
1.5-2.0, which is adequate for our purposes. 

To estimate the emergent luminosity from the disk, wc assume that the central star is the only 
radiation source. We follow Kenyon & Hartmann (1987) and assume a spherical, limb-darkened 
star with radius Rq, luminosity Lq, and limb-darkening coefficient eo = 0.6. For a point P at the 
outer boundary of annulus i with height hp above the disk midplane, rays from the star enter the 
annulus at a scale height /ij„ above (below) the midplane. We compute the length I of the path 
through the disk and derive the optical depth along this path as Tp = {l/Aai)Ti, where Aoj is 
the width of the annulus. The radiation absorbed along this path is e~'^p/o, where Iq is the flux 
incident on the boundary of the annulus. Numerical integrations over the stellar surface and the 
vertical extent of an annulus yield the amount of flux absorbed by each annulus, which we convert 
to relative surface brightness. A final numerical integration over the radial extent of the disk yields 
the ratio of the disk luminosity to the stellar luminosity, Ld/Lq. 

Although this algorithm is an improvement over our previous optically thin treatment of radia- 
tive transfer through a planetesimal disk, it has several limitations. The optical depth calculation 
does not distinguish between absorption and scattering. For a grain albedo ui, the luminosity in 
scattered light is roughly ujLd/Lq; the thermal luminosity is roughly {1 — lo)Ld/Lq. The algorithm 
also ignores multiple scattering and absorption. Because the vertical optical depth of the disk is 
~ 10~^ or smaller in most cases, this approximation is satisfactory. 
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To download the figures for this paper, including two color stills, go to 
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To view the two animations described in the paper, go to 
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